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Resumen 


El conocimiento global de las consecuencias del cambio de vegetación y 
uso de suelo en las fases del ciclo hidrológico y la calidad del suelo aún 
está en ciernes. Esta tarea es particularmente desafiante en las regiones 
tropicales donde las tasas de deforestación alcanzan los niveles máximos 
reportados y el monitoreo es limitado. Tal es el caso de México, donde el 
cambio de cobertura y la erosión hídrica afectan a más de la mitad del 
territorio nacional. En este contexto, se usó el modelo SWAT para estudiar 
la erosión hídrica asociada con las variaciones en el caudal condicionadas 
por la trayectoria de cambio de la cobertura vegetal en seis periodos entre 


1986 y 2018 en dos subcuencas del río Papaloapan. Basados en la 
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evaluación combinada de tres estadísticos de eficiencia y dos de 
incertidumbre, 6 de los 10 periodos calibrados mostraron buenos ajustes 
en la simulación de la descarga con KGE > 0.70 y P-factor > 0.70. Los 
resultados mostraron a mediano plazo una disminución en la producción 
de sedimentos en las dos subcuencas estudiadas, posiblemente 
relacionada con la recuperación de los bosques en las cabeceras de las 
cuencas. A pesar de los buenos resultados de las simulaciones, el análisis 
detallado de las fuentes de incertidumbre en regiones con monitoreo 
limitado es tan importante como la validación con datos medidos de buena 
calidad para mejorar el desempeño y la confianza de los modelos y 


soportar de forma adecuada la toma de decisiones. 


Palabras clave: erosión, caudal, SWAT, deforestación, cambio de uso de 


suelo, trópicos, eficiencia, incertidumbre. 


Abstract 


Global knowledge on the consequences of changes in vegetation and land 
use in the phases of the hydrological cycle and soil quality ¡is still incipient. 
This task ¡is particularly challenging in tropical regions, where 
deforestation rates reach their highest reported levels and monitoring is 
limited. Such is the case of Mexico, where changes in coverage and water 
erosion affect more than half of the national territory. In this context, the 


SWAT model was used to study water erosion associated with variations 
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in flow conditioned by the trajectory of change in vegetation cover in six 
periods, between 1986 and 2018, in two sub-basins of the Papaloapan 
River. Based on the combined evaluation of three efficiency and two 
uncertainty statistics, six of the ten calibrated periods showed good fitting 
in discharge simulation with KGE > 0.70 and P-factor > 0.70. In the 
medium term, the results showed a decrease in sediment yield in the two 
study sub-basins, possibly related to the recovery of the forests in the 
basins' headwaters. In spite of the good results of the simulations, the 
detailed analysis of the sources of uncertainty in regions with limited 
monitoring is as important as the validation with good-quality measured 
data, in order to improve the performance and reliability of the models 


and to support adequate decision-making. 


Keywords: Erosion, flow, SWAT, deforestation, land use change, tropics, 


efficiency, uncertainty. 
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Introducción 


La erosión hídrica es un proceso de la superficie terrestre de origen 
natural y antrópico cuya magnitud a escala mundial se estima en un rango 
de 20-30 gigatoneladas anuales (FAO € ITPS, 2015). En las regiones 
tropicales del sur global ocurren las mayores tasas de deforestación 
(Borrelli et a/l., 2017) promovidas por el avance de la frontera 
agropecuaria, la extracción de recursos o la urbanización (Montes-León, 
Uribe-Alcántara, 8 García-Celis, 2011). Durante la década de 2000 a 
2010, las transiciones de coberturas de bosques a otros usos de suelo 
avanzaron a un ritmo de 0.46 % en América Latina y el Caribe (FAO, 
2011), y a 0.3 % en México (Blackman, Goff, 8 Rivera-Planter, 2018); 
ambas tasas superan el promedio mundial de 0.13 % (FAO, 2011). 
Además de la actividad humana, la erosión hídrica en los trópicos se 
agudiza debido a los regímenes de lluvias estacionales de gran intensidad 
(Rangel, Jorge, Guerra, 8: Fullen, 2019). De esta manera, más del 60 % 
de los suelos de México presenta algún nivel de erosión por degradación 
hídrica, la cual se agudiza a una tasa anual de 0.5 % en promedio (Borrelli 
et al., 2017). Estos graves procesos de erosión hídrica se agudizan 


especialmente en las etapas subsecuentes a los cambios de cobertura y 
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uso de suelo (Cardoza-Vázquez et al., 2007; Bolaños-González et al., 
2016). 


Los servicios ambientales, el ciclo hidrológico y el suelo son 
afectados por los cambios en la cobertura a múltiples escalas espaciales 
y temporales dependiendo de la fisiografía de la cuenca, de la variabilidad 
climática y de la trayectoria de uso del territorio estudiado (Martinez, 
Guyot, Filizola, 8 Sondag, 2009; Arheimer 8 Lindstróm, 2019). Por 
ejemplo, la vegetación y el uso de suelo son factores importantes en el 
comportamiento hidrológico de las cuencas, pues determinan los procesos 
de generación del escurrimiento superficial y subsuperficial (Arheimer € 
Lindstróm, 2019; Couto et a/., 2019; Liu, Zhang, Kou, 4 Zhou, 2017). A 
escala regional, se incrementa el volumen de la escorrentía y de los 
sedimentos transportados, con alteraciones morfológicas a la red fluvial y 
el régimen hidrológico de las cuencas, y aumenta el riesgo de 
inundaciones en planicies (Chotpantarat 8 Boonkaewwan, 2018; Wang, 
Ge, Yu, Wang, 84 Xu, 2017; Yang et al., 2018). 


En los países tropicales, el monitoreo hidrológico se realiza con 
redes de medición mínimas (Vórósmarty, Lévéque, 8: Revenga, 2005), lo 
que se traduce en registros históricos limitados (Carvalho, 2019; Delmas, 
Cerdan, Cheviron, € Mouchel, 2011; Warrick, 2014). La convergencia 
entre la escasez de datos históricos, la necesidad de información para la 
toma de decisiones, el avance en las tecnologías de cómputo y la 


producción de conocimiento científico favorecen la aplicación de modelos 
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hidrológicos como las herramientas más utilizadas de investigación, 
planificación y manejo ambiental (Fakhri, Dkohaki, Eslamian, Fazeli- 
Farzani, € Reza-Farzaneh, 2014; Jajarmizadeh, Harun, 8 Salarpour, 
2012; Jorgensen 8 Fath, 2011). De igual manera, los sistemas de 
información geográfica han aumentado las capacidades para gestionar 
grandes cantidades de información, con ventajas como el despliegue de 
los resultados distribuidos espacialmente y la generación de escenarios 
(Alatorre 8. Beguería, 2009; Goodchild, 1996). 


La herramienta de evaluación de suelo y agua (SWAT, por sus siglas 
en inglés) es un modelo conceptual con una base física que estima la 
producción de agua, sedimentos, nutrientes, contaminantes y patógenos 
a escala de cuenca a partir de datos espacialmente distribuidos acerca de 
la topografía, el suelo, la vegetación y el uso de suelo, las prácticas de 
manejo y el clima (Douglas-Mankin, Srinivasan, € Arnold, 2010). El 
funcionamiento de este modelo está definido bajo los principios de masa 
y de momentum: el primero se relaciona con el balance hídrico, de 
sedimentos, nutrientes o contaminantes y el segundo se enfoca sobre 
todo en el desplazamiento de estos flujos a través de la red de drenaje 
de la cuenca. Entre los principales ámbitos de aplicación del modelo SWAT 
destacan aspectos de la calidad de agua como la contaminación difusa, la 
producción de sedimentos, la simulación de la escorrentía y los impactos 


hidrológicos del cambio ambiental (Wang et a/., 2019). 
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En el modelo SWAT, la dinámica de las variables de interés es 
modelada a partir de datos de entrada fisiográficos, climáticos y 
antrópicos espacialmente distribuidos, de manera que es de utilidad en 
cuencas no aforadas o en la evaluación de escenarios a largo plazo en 
función de variaciones de alguna de las variables de entrada (Neitsch, 
Arnold, Kiniry, 8 Williams, 2011); en este caso de estudio será en función 
de los cambios de cobertura y uso de suelo. El empleo del modelo SWAT 
en países con monitoreo limitado es extensa. En distintas regiones de 
África oriental, Guzha, Rufino, Okoth, Jacobs y Nóbrega (2018) 
identificaron 15 estudios que utilizaron SWAT para evaluar el impacto del 
uso de suelo en la escorrentía superficial, la descarga y el flujo base. Por 
su parte, Tan, Gassman, Srinivasan, Arnold y Yang (2019) revisaron 126 
casos que estudiaron diferentes temas de cambio de uso de suelo, cambio 
climático, prácticas de manejo o calidad del agua con SWAT en el sureste 
asiático. 

En la cuenca del río Papaloapan, Bello, Gómez, Magaña, Graizbord 
y Rodríguez (2009), y Ruiz-Fernández et al. (2014) señalaron que los 
problemas de azolvamiento del sistema de lagunas y humedales de la 
cuenca baja están asociados con una gran cantidad de sedimentos 
exportada por los afluentes aguas arriba durante los últimos 40 años. Por 
su parte, Pérez-Vega y Ortiz-Pérez (2002) relacionaron positivamente el 


incremento en el escurrimiento de las subcuencas de los ríos San Juan y 
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Tesechoacán con el incremento de la deforestación en la cuenca del río 


Papaloapan entre 1973 y 1993. 


Con la perspectiva de que los procesos ambientales son históricos 
—en cuanto producto de un legado histórico, pero también resultado de 
una conjugación de variables en el presente (Reid 1998; Verstraeten 
Lang, 8 Houben, 2009)— y susceptibles de estudio con los recursos 
disponibles, se usó el modelo SWAT para reconstruir la trayectoria de la 
erosión potencial asociada con el cambio de uso de suelo en seguimiento 


a los objetivos específicos: 


Es Estimar las variaciones de la erosión hídrica y el caudal con 
respecto al cambio de uso y cobertura de suelo de 1985 a 2018 
en dos subcuencas del río Papaloapan. 

2. Estimar la distribución espacial del potencial erosivo en el área 


de estudio. 


Área de estudio y bases de datos 
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Las subcuencas de los ríos San Juan (SIN) y Tesechoacán (TES) 
constituyen el 26 % de los 46 720 km? de área total de la cuenca del río 
Papaloapan que se localiza en la vertiente del Golfo de México entre los 
190 N 970 W y los 170 N 940 W (Figura 1). 
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Figura 1. Localización de las subcuencas estudiadas (derecha superior) 


con elevación (izquierda) y malla de datos del Sistema Global de 
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Asimilación de Datos Terrestres (GLDAS, por su acrónimo en inglés) 


(derecha inferior). 


Alrededor del 55 % del área estudiada tiene pendientes inferiores a 
15 % y el 45 % del área restante pendientes mayores a 15 % (FAO, 
2009). Siguiendo la clasificación de Kóppen modificada para México por 
García (2004), en la región se identifican tres grupos climáticos: (A) 
tropical lluvioso, en las planicies próximas al nivel del mar; (B) climas 
secos con régimen de lluvias en verano y (C) templado húmedo en las 
áreas de las cuencas altas (hasta 3 400 msnm). Se identificaron 12 tipos 
de vegetación y uso de suelo agrupados en cinco clases: bosques 
templados, selvas tropicales, matorrales y chaparrales, pastos y 
agricultura (Alavez-Vargas, Birkel, Corona, 8 Breña-Naranjo, 2021) 
(Figura 2). 
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Figura 2. Vegetación y uso de suelo de 1986 a 2018, representación 


geográfica y porcentajes por subcuenca estudiada. 
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En este estudio se usó el modelo digital de elevación de la Misión 
Topográfica de Radar Aerotransbordado (SRTM, por su acrónimo en 
inglés) con una resolución espacial de 90 m (Jarvis, Reuter, Nelson, 8 
Guevara, 2008). Esta información ráster fue reclasificada en rangos de 
pendientes en función de la geomorfología de la cuenca y la clasificación 
de la FAO (2009) para formar dos clases: 0-15 % y mayor a 15 %. Los 
conjuntos vectoriales correspondientes a la red de drenaje (INEGI, 2010) 
provienen de la información geoespacial temática producida por el 
Instituto Nacional de Estadística, Geografía e Informática (INEGI) de 
México. La información edafológica requerida por el modelo fue derivada 
de los conjuntos vectoriales relativos a las unidades edafológicas (INEGI, 
2007) y los perfiles edafológicos (INEGI, 2013) además de dos softwares 
para estimar otras características no disponibles en las fuentes 
consultadas. Con Soil Water Characteristics (Saxton € Rawls, 2006) se 
estimaron la conductividad hidráulica saturada, la densidad aparente y la 
infiltración básica, mientras que con NumCur 10 IE se identificaron los 
grupos hidrológicos de suelo. Los tipos de vegetación y uso de suelo se 
clasificaron siguiendo el procedimiento descrito en Alavez-Vargas et al. 
(2021) a partir de imágenes Landsat descargadas de Google Earth Pro 
(Wuthrich, 2006) procesadas en QGIS v.3.4.5 Madeira (QGIS 
Development Team, 2019) para los años 1986, 1993, 1998, 2003, 2011 
y 2018 (Figura 2) y posteriormente reclasificadas conforme a la 
nomenclatura usada por el modelo. Los datos de temperatura, 


precipitación, velocidad del viento y radiación solar diarios (Li, Beaudoing, 
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8: Rodell, 2018) y trihorarios (Beaudoing € Rodell, 2015; Beaudoing 8 
Rodell, 2016) con resolución espacial aproximada de 25 km? para el 
periodo de 1980-2018 se obtuvieron del portal del Sistema Global de 


Asimilación de Datos Terrestres (GLDAS, por su acrónimo en inglés) 


(Rodell et a/., 2004), en una malla de 63 puntos sobre el área estudiada 


(Figura 1). Finalmente, las mediciones de caudal diario de 1985 a 2014 


se extrajeron del Banco Nacional de Datos de Aguas Superficiales 


(BANDAS por sus siglas en español) (Conagua, 2017). Los productos, 


resoluciones y fuentes consultadas para obtener los datos de las 10 


variables biofísicas y climáticas recabadas se resumen en la Tabla 1. 


Tabla 1. Caracterización de los datos de entrada. 


Producto y resolución 


Variable Resolución Fuente 
temporal 
Modelo digital de : : 
ne SRTM, formato ráster 90 m Jarvis et al. (2008) 
elevación 
: Red hidrográfica, formato Escala 
Red hidrica INEGI (2010) 
vectorial 1:50 000 
Unidades edafológicas INEGI (2007) 
, ! Escala 
Edafología Carta de perfiles 
1:250 000 


Open Access bajo 


edafológicos 


la licencia CC BY-NC-SA 4.0 
(https: //creativecommons.org/licenses/by-nc-sa/4.0/) 


INEGI (2013) 
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- Producto y resolución E 
Variable Resolución Fuente 

temporal 
Escala 
1:250 000 
Imágenes Landsat, 
o formato ráster, años 

Vegetación 30m Wuthrich (2006) 


Temperatura 


1986, 1993, 1998, 2003, 


2011, 2018 


GLDAS Noah Land 
Surface Model L4 V 2.0 
Periodo: 1980-01-01 a 

2010-12-31 


GLDAS Noah Land 


Tres horas, 
0.259 x 0.250 


Tres horas, 


Beaudoing y Rodell 
(2015) 


Beaudoing y Rodell 


Surface Model L4 V 2.1 0.259 x 0,250 
(2016) 
Periodo: 2011-01-01 a 
2018-12-31 
Precipitación GLDAS Catchment Diario Li et al. (2018) 
Velocidad del surface Model L4 V 2.0 0.259 x 0.250 
viento Periodo: 1980-01-01 a 


Open Access bajo la licencia CC BY-NC-SA 4.0 
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Producto y resolución 


Variable Resolución Fuente 
temporal 
2010-12-31 
3 horas Beaudoing y Rodell 
GLDAS Noah Land 0.259 x 0,259 (2016) 


Surface Model L4 V 2.1 
Periodo: 2010-01-01 a 
2018-12-31 


Radiación solar 


Localización de las 


estaciones hidrométricas 
Estaciones A , 
] 0 Garro (rio Tesechoacán) y Conagua (2017) 
hidrométricas : 
Cuatotolapan (rio San 


Juan) 


Registros diarios del 


BANDAS 
Caudal medido diarios Conagua (2017) 
Periodo 1985-01-01 a 


2014-12-31 
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Métodos 


Configuración del modelo SWAT 


La modelación en SWAT consta de cinco etapas (Figura 3). Al inicio, el 
modelo digital de elevación y la red fluvial guían la delimitación de las 
subcuencas de interés. Después se incorporan las variables con las que el 
modelo define las unidades de respuesta hidrológica (HRUs): pendiente, 
tipo de suelo, tipo de vegetación y uso de suelo. Estas dos primeras 
etapas, por lo tanto, se enfocan en la delimitación física de las unidades 
(subcuencas) y subunidades de análisis (HRUS) a través de las cuales, en 
la etapa siguiente, se agrupa la información de las variables de entrada y 
se calculan las respuestas de los componentes del balance hídrico 
(Ecuación (1)) y las variables de interés (la erosión en este caso, Ecuación 
(2)) (Neitsch et a/., 2011). 
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INPUTS FASE DE LA MODELACIÓN 


Modelo digital 


de elevación 


Red hídrica 


Pendiente 


Delimitación de subcuencas 


Características 
edafológicas 


Unidades hidrológicas 


Vegetación y de respuesta 
uso de suelo 


Temperatura 


Precipitación 


Velocidad Calibración 
del viento 
Parametrización 


Radiación solar 
Caudal medido Escenario e 
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Figura 3. Flujo de trabajo de la modelación en SWAT: preparación de 


datos de entrada, aplicación y calibración del modelo. 


En la tercera etapa se introdujeron al modelo las variables 
climáticas: precipitación, temperaturas máximas y mínimas, radiación 
solar y velocidad del viento. Estas variables representan en conjunto las 
cantidades de agua y energía que gobiernan el balance hídrico. En 
concreto, estos datos son usados por el modelo para resolver la ecuación 
de balance hídrico (Ecuación (1)) a escala de HRUs por subcuenca 
(Neitsch et al., 2011): 


SW, = SWo + Yi; Lan = Oief — Ea — Wseep — 0) (1) 


Donde SW: es el contenido final de agua en el suelo (mm); SWo, el 
contenido inicial de agua en el suelo en el día ¡ (mm); t, el tiempo en días; 
Raay, la precipitación en el día ¡ (mm); Qsurr, la escorrentía superficial en 
el día ¡ (mm); Ex, la evapotranspiración en el día ¡ (mm); Wseep, la 


percolación en el día ¡ (mm), y Qgw es el flujo de retorno en el día ¡ (mm). 


A su vez, la erosión en HRUs por subcuenca es calculada con la 
ecuación universal de la pérdida de suelo modificada (MUSLE, Ecuación 
(2)) (Williams, 2012): 
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0.56 
sed = 11.8 (Osirg X peak % ATC pru) X Kysie X CusiE X Pusie X LSysig X CFRG 


(2) 


Donde sed es el sedimento producido diario (ton); Qsurr, la lámina 
de escorrentía superficial diaria (mm ha); qpeax, el caudal pico (m* s”1); 
areahru, el área de la HRU (ha); Kuse, el factor de erodabilidad del suelo; 
Cus.e, el factor de cobertura vegetal y manejo; Puse, el factor de prácticas 
de soporte; LSus.e, el factor topográfico, y CFRG es el factor de 


fragmentación del suelo. 


El resultado de la tercera etapa es la simulación inicial del escenario 
en turno. Esta simulación es calibrada con datos observados de las 


variables analizadas, como se describe en la sección siguiente. 


Simulación inicial y calibración 


Este estudio se realizó en QGIS versión 2.6.1 Brighton (QGIS 


Development Team, 2014) con la extensión QSWAT versión 1.8 (Dile, 
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Srinivasan, 8 George, 2019). Los datos de entrada fueron preparados en 
el formato requerido por el software para el cálculo de las diferentes 
variables involucradas en las ecuaciones (1) y (2). Se usaron las 
ecuaciones establecidas por defecto por el modelo (ver Arnold et al., 
2011), excepto la evapotranspiración, que se calculó con el método de 
Hargreaves, Hargreaves y Riley (1985). Los valores de los parámetros 
físicos (es decir, la elevación, la red de drenaje, la pendiente y las 
características edafológicas) fueron constantes para los seis periodos 
evaluados en cada subcuenca. Las coberturas de vegetación y uso de 
suelo fueron específicas para cada escenario. Por su parte, la información 
climática introducida tuvo resolución diaria. Debido a la ausencia de bases 
de datos sobre las prácticas de conservación de suelo en el área de 
estudio, se consideró el factor de prácticas de soporte Pus. igual a 1, es 


decir, sin prácticas de conservación. 


Para el periodo de 1986 a 2018 se realizaron seis modelaciones por 
subcuenca correspondientes a los periodos definidos por los escenarios 
de vegetación y uso de suelo, según se muestra en la Tabla 2. La 
irregularidad de los intervalos evaluados se debe a las decisiones 
metodológicas realizadas para asegurar el mosaico de imágenes de la 
mejor calidad posible para caracterizar la cobertura y el uso de suelo en 


cada escenario. 
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Tabla 2. Periodos de calentamiento, modelación y calibración en las tres 


subcuencas estudiadas. 


Escenario de Calibración 
vegetación y uso | Modelación | Calentamiento de caudal 
de suelo mensual 
1986 1982-1986 1982-1983 1985-1986 

1993 1985-1993 1985-1986 1987-1993 

1998 1992-1998 1992-1993 1994-1998 

2003 1997-2003 1997-1998 1999-2003 

2011 2002-2011 2002-2003 2004-2011 

2018 2010-2018 2010-2011 2012-2014 


Para cada modelación se definió un periodo de calentamiento de dos 
años como sugiere Arnold et a/. (2012). La escala temporal de la 


modelación fue mensual. 


La calibración del caudal fue realizada con 30 años de datos 
observados (1985-2014) en las estaciones hidrométricas de las 
subcuencas San Juan y Tesechoacán. La erosión simulada no fue calibrada 
debido a la ausencia de datos de erosión medida en las dos subcuencas 


estudiadas. 
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Los análisis de sensibilidad, incertidumbre y calibración de los 
modelos se hicieron siguiendo las recomendaciones de Abbaspour, 
Vaghefi y Srinivasan (2017) con el algoritmo SUFI-2 disponible en el 
software SWAT-CUP (Abbaspour, 2015). El algoritmo SUFI-2 se enfoca en 
la búsqueda del grupo de parámetros con mejor ajuste y menor 
incertidumbre respecto a los datos observados mediante el muestreo por 
el método del hipercubo latino (Beven, 2012), con el objetivo de cumplir 
dos condiciones: que más del 90 % de los datos se localicen en el área 
de predicción del 95 % (95PPU, por el acrónimo en inglés de 95 % 
prediction uncertainty) y que la proporción en la diferencia promedio entre 
los valores límite del 95PPU y la desviación estándar sea menor a 1 


(Abbaspour, Johnson, € van Genuchten, 2004). 


Para cada escenario y subcuenca se realizaron al menos tres 
iteraciones de 600 a 1 000 simulaciones. Se utilizó el criterio de Kling 
Gupta (KGE) como función objetivo y además se emplearon los 
estadísticos de eficiencia de Nash-Sutcliffe (NSE), el porcentaje de sesgo 


(PBIAS) y el coeficiente de determinación (R2). 


El criterio de Kling Gupta evalúa la bondad de ajuste del modelo 
mediante el cálculo de las distancias euclidianas entre la correlación, el 
sesgo y la variabilidad del valor simulado respecto al valor observado 
(Ecuación (3) y Ecuación (4)) como una alternativa al criterio de Nash- 
Sutcliffe que toma como referencia la diferencia de 1 menos el cociente 


del error cuadrático medio y la varianza (Ecuación (5)) para mostrar la 
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proporción de la variabilidad de las observaciones explicada por la 


simulación (Gupta, Kling, Yilmaz, € Martínez, 2009): 


KGE =1-ED (3) 


ED == D7+(a-D7+ (6-17 (4) 


Donde ED es la distancia euclideana desde el punto ideal; r, el 
coeficiente de correlación lineal entre datos simulados y observados; a,la 
variabilidad dada por el cociente entre las desviaciones estándar de los 
datos simulados y los observados, y finalmente f es el sesgo dado por el 


cociente entre los flujos medios simulado y observado: 


E E (ros) pe MSE 
NSE = 1- === =1-—3 (5) 
¿=1(%o,—Mo) a] 


Donde n es el número total de pasos de tiempo; Xs+, el valor 
simulado al tiempo t, y o y coson la media y la desviación estándar de 


los valores observados. 


El porcentaje de sesgo evalúa la tendencia media de los datos 


simulados a sobreestimar o subestimar respecto a los observados. Los 
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factores P y R son estadísticos de incertidumbre respecto al intervalo de 
predicción 95PPU, donde el P-factor se refiere al porcentaje de datos 
observados situados en la banda de predicción 95PPU y el R-factor indica 
la amplitud de la banda de predicción 95PPU (Abbaspour, 2015), que 


cuanto más estrecha indica una menor incertidumbre (Abbaspour, 2005). 


A partir de la mejor simulación inicial y los análisis de sensibilidad, 
el software sugiere una serie de parámetros significativos para realizar la 
calibración (Tabla 3). En cada una de las subcuencas analizadas se calibró 
el caudal para cada uno de los seis escenarios, con un número variable 
entre 4 y 8 parámetros (Tabla 4 y Tabla 5), dependiendo de la sensibilidad 


de cada simulación a las variaciones en los parámetros involucrados. 
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Tabla 3. Parámetros usados en la calibración de la descarga de las tres 


subcuencas. 


Acrónimo 


R_CN2.mgt 


Descripción del parámetro 


Número de curva para la condición de 
humedad II 


Componente 
hidrológico 


Escorrentía 
superficial 


V_ALPHA_BF.gw 


Factor alfa de la curva de recesión del 
hidrograma (1/días) 


V_GW_DELAY.gw 


Tiempo de retraso del agua subterránea (días) 


V_GWQMN.gw 


Umbral de profundidad del agua en el acuífero 
somero requerido para generar el flujo base 
(mm H20) 


Agua subterránea 


R_SOL_AWC(..).sol 


Capacidad de agua disponible en el suelo (mm 
H20/mm de suelo) 


Agua del suelo 


R_ESCO.hru 


Factor de compensación de evaporación del 
suelo a escala de HRU 


R_ESCO.bsn 


Factor de compensación de evaporación del 
suelo a nivel de cuenca 


V_GW_REVAP.gw 


Coeficiente de aguas del acuífero somero que 
regresan a la zona radicular 


V_REVAPMN.gw 


Umbral de profundidad del agua en el acuífero 
somero que regresa a la zona radicular (mm 
H20) 


Evapotranspiración 
potencial y real 
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Tabla 4. Rangos iniciales y finales de los parámetros usados en la 


OPEN a ACCESS 


1) Check for updates 


calibración del caudal y estadísticos de eficiencia de la cuenca TES para 


cada uno de los periodos modelados. 


Rango inicial (mín/máx) 


Núm. Acrónimo Periodo 

1985-1986 1987-1993 1994-1998 1999-2003 2004-2011 2012-2018 
1 R_CN2.mgt -0.1/ 0.11 -0.16/ -0.12 -0.2/ -0.07 -0.2/ -0.12 -0.2/ -0.15 -0.2/ -0.18 
2 V_ALPHA_BF.gw 0.23/ 0.36 0.26/ 0.35 -0.01/ 0.95 -0.01/ 1.02 -0.05/ 0.2 0.29/ 0.4 
3 V_GW_DELAY.gw | 98.03/ 192.86 58.75/ 95.52 -3.08/ 262.89 | 277.78/ 541.58 | 20.21/ 152.76 | 123.94/ 203.7 
4 V_GWQMN.gw 1.95/ 2.42 111.32/ 178.11 -1.38/ 0.34 -9.17/ 144.7 1.15/ 2.11 0.53/ 1.2 
5 R_SOL_AWC(..).sol -0.34/ 0.72 0/ 0.06 -0.34/ -0.14 0/ 0.77 0.02/ 0.48 -0.88/ -0.41 
6 R_ESCO.hru -0.12/ 0.51 
7 R_ESCO.bsn 0.31/ 0.72 -0.41/ -0.14 0.27/ 1.61 0.63/ 0.94 -0.02/ 0.26 

Ver en la Tabla 3 la descripción del parámetro representado con el acrónimo. 
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Tabla 5. Rangos iniciales y finales de los parámetros usados en la 


calibración del caudal y estadísticos de eficiencia de la cuenca SJN para 


cada uno de los periodos modelados. 


Rango inicial (mín/máx) 


Núm Acrónimo Periodo 

1985-1986 1987-1993 1994-1998 1999-2003 2004-2011 2012-2018 
1 R_CN2.mgt -0.09/ 0.15 -0.09/ 0 -0.2/ 0.05 -0.2/ 0.06 -0.19/ -0.16 -0.2/ 0.04 
2 V_ALPHA_BF.gw 0.32/ 1.03 0.02/ 0.39 -0.11/ 0.66 -0.4/ 0.55 0.39/ 0.61 0.45/ 1.4 
3 V_GW_DELAY.gw | 207.91/ 587.09 | -8.05/ 122.06 | -139.63/ 262.63 | 49.77/ 325.23 | 290.53/ 355.01 | -107.24/ 272.24 
4 V_GWQMN.gw 0.1/ 1.4 1.48/ 2.15 0.94/ 2.96 44.59/ 140.41 | 121.19/ 156.32 -0.65/ 1.15 
5 R_SOL_AWC(..).sol -0.57/ 0.33 0.51/ 0.66 0.08/ 0.71 
6 R_ESCO.bsn -1.3/ -0.45 0.45/ 0.67 0.25/ 0.8 
7 V_GW_REVAP.gw 0.16/ 0.21 
8 V_REVAPMN.gw 17.94/ 159.06 

Ver en la Tabla 3 la descripción del parámetro representado con el acrónimo. 
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Resultados 


El cambio de uso de suelo 


La superficie de la subcuenca TES cubierta por bosques templados 
aumentó de 16 % en 1986 a 26 % en 2003 y luego declinó nuevamente 
en 2018 (Figura 2, gráfica de áreas TES). Por su parte, las selvas 
mostraron una pérdida de 50 % de su extensión de 1986 a 2018, con un 
marcado declive en 2003 y posteriormente una recuperación en los 
siguientes 15 años. Las áreas ocupadas por chaparral y matorral rosetófilo 


fueron inferiores al 0.25 % en todos los escenarios. 


Las variaciones de las coberturas vegetales naturales están 
relacionadas directamente con variaciones en las prácticas productivas. 
La recuperación de los bosques templados se debió sobre todo al 
abandono de la agricultura de temporal en las partes altas entre 1986 y 
2003, mientras que el declive en el periodo de 2003 a 2018 indica 
prácticas humanas que propician la deforestación. A su vez, la reducción 


de la extensión de la selva en las partes bajas ocurrió por el desmonte 
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realizado para la expansión de las áreas de pastizal (con fines pecuarios), 
que llegaron a ocupar la mitad de la subcuenca en 2003. En los años 
subsecuentes, el abandono de los pastizales redujo su extensión hasta 
ocupar solo el 30 % de la cuenca, permitiendo la recuperación de la 
vegetación secundaria de las selvas en esos espacios. De 1986 a 2003, 
las áreas agrícolas se redujeron de modo paulatino de 12 % a menos del 
1 % del área; luego se ampliaron hasta ocupar el 20 % de la subcuenca 
en 2018. 


La subcuenca San Juan es la más meridional y de menor altitud (< 
1 000 msnm) de manera que las áreas cubiertas con bosques templados 
son menores (< 5 %, excepto en 2003 y 2011) que en TES (Figura 2, 
gráfica de áreas SJN). Por el contrario, las selvas húmedas llegaron a 
cubrir el 55 % del área en 1993, seguida por un proceso de degradación 
que redujo esta extensión a una cuarta parte en 2003 y una fase de 
recuperación en las últimas dos décadas. En todos los escenarios, el 


chaparral representó menos del 0.05 % del área (> 400 ha). 


De manera semejante a TES, el área perdida por bosques y selvas 
de 1986 a 2003 fue ganada por los pastizales que llegaron a cubrir el 68 
% de la subcuenca en 2003 y reducida después a 39 %. Por su parte, las 
áreas agrícolas han incrementado el área ocupada de 10 % en 1986 a 25 
% en 2018, con el punto más bajo alcanzado en 2003. Sin embargo, como 
se discute más adelante en relación con la calidad de los datos de entrada 


(ver sección Las fuentes de incertidumbre y los alcances de las 
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modelaciones), los cambios drásticos registrados en 2003 respecto al 
escenario previo (1998) y posterior (2011) sugieren la necesidad de 
verificar los resultados obtenidos para ese año con un set de imágenes 


diferente. 


La modelación hidrológica 


Los resultados que se presentan en esta sección y en la siguiente son 
producto de la parametrización de los modelos con el mejor ajuste 
obtenido según los estadísticos de eficiencia. En la Tabla 4 (TES) y la 
Tabla 5 (SJN) se muestran los rangos de incertidumbre iniciales dados por 
defecto por el modelo en referencia a los cuales se ajustaron los valores 


finales de cada parámetro para cada subcuenca. 


Las gráficas con los datos de caudal parametrizados con el mejor 
ajuste y los intervalos de predicción para cada periodo se presentan en la 
Figura 4 y la Figura 5, junto con los respectivos estadísticos de eficiencia 


e incertidumbre. 
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Figura 4. Modelación del caudal de 1986-2014 (calibrado) y de 2015- 


2018 (simulado) en la estación hidrométrica “Garro” de la subcuenca 
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Tesechoacán (TES). Cada gráfica (inciso) corresponde a un periodo 


modelado (ver Tabla 2). 


a) 1985-1986 b) 1987-1993 c) 1994-1998 
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Figura 5. Modelación del caudal de 1986-2014 (calibrado) y 2015-2018 
(simulado) en la estación hidrométrica “Cuatotolapan” de la subcuenca 
San Juan (SIN). Cada gráfica (inciso) corresponde a un periodo 
modelado (ver Tabla 2). 


Independientemente de los cambios de vegetación y uso de suelo 
entre los diferentes escenarios por subcuenca, los mismos seis 


parámetros presentaron el mejor ajuste en TES (Tabla 4). 


El porcentaje de sesgo indica que los resultados son buenos cuando 
se encuentran en el rango de + 15 % (Moriasi et al., 2007). Abbaspour 
(2015) sugiere que, en la calibración de descargas, los valores de P-factor 
> 70 % se consideran óptimos. El mismo autor señala que el valor de R- 
factor debe ser menor a 1.5 y para esta subcuenca todos los periodos 
presentaron valores entre 0.31 y 0.58. 


TES presenta valores de eficiencia en todos los periodos calibrados 
entre 0.64 y 0.90, excepto en el periodo 1999-2003 (KGE = 0.58, NSE = 
0.14, R? = 0.35) (Figura 4). El bajo desempeño del periodo 1999-2003 se 
verifica en el bajo valor P-factor = 0.28, el coeficiente de sesgo 
comparativamente alto (PBIAS = 7.5), así como de manera visual (Figura 
4 inciso d). Los estadísticos de incertidumbre confirman el ajuste 
satisfactorio en los periodos 1994-1998, 2004-2011 y 2012-2018 (valores 
de P-factor mayores a 70 % en los tres casos) e indican que el periodo 


1985-1986 queda ligeramente por debajo del umbral de buen desempeño 
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(P-factor = 0.63). El primero y tercer periodos sobreestiman ligeramente; 
el resto de los periodos presenta una ligera subestimación. A pesar de 
que el periodo de los datos simulados de 1987-1993 presentó los mejores 
valores de KGE, NSE y R? (> 0.80), apenas el 40 % de los datos 
observados en 1987-1993 coinciden con el intervalo 95PPU (Figura 4 


inciso b). 


En SIN, la escorrentía superficial y el agua subterránea, con 1 
(R_CN2.mgt) y 3 (V_ALPHA_BF.gw, V_GW_DELAY.gw, V_GWQMN.gw) 
parámetros, respectivamente, explicaron el comportamiento del 
rendimiento hídrico de 1986 a 1998; mientras que de 2003 a 2018, el 
agua del suelo (R_SOL_AWC(..).sol) y la evapotranspiración potencial y 


real tuvieron mayor relevancia (Tabla 5). 


De manera semejante a TES, en SJN (Figura 5), el periodo 1999- 
2003 presenta los valores de eficiencia más bajos (< 0.54). Por el 
contrario, los periodos 1987-1993, 1994-1998 y 2012-2018 mostraron 
los mejores valores de desempeño tanto en los estadísticos KGE (0.90, 
0.78, 0.89, respectivamente), NSE (0.81, 0.57 y 0.83, respectivamente) 
y R? (0.83, 0.62 y 0.85, respectivamente) como en los factores de 
incertidumbre P-factor (0.82, 0.72, 0.75, respectivamente) y R-factor 
(0.72, 0.72, 0.75, respectivamente). Si bien los periodos 1985-1986 y 
2004-2011 presentan buenos ajustes, el factor P indica un nivel de error 


mayor con valores de 0.54 y 0.57, respectivamente. 
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La modelación de la erosión potencial 


Los valores de erosión potencial estimados en esta sección corresponden 
a los promedios multianuales de los seis periodos de caudal modelados. 
La erosión hídrica media por escenario por subcuenca oscila en el rango 
de 69 a 158 ton ha? año"? en la subcuenca Tesechoacán, y entre 57 y 
162 ton ha"! año"? en San Juan (Tabla 6). De acuerdo con estas cifras, el 
volumen erosionado en Tesechoacán ha disminuido a 44 % del volumen 


estimado en 1986 y en San Juan a 38 % de la estimación inicial. 


Tabla 6. Erosión hídrica promedio multianual (ton ha"? año"1) por 


periodo por subcuenca. 
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Erosión hídrica media multianual (ton ha! año"*) 
Subcuenca Periodo 
1985-1986 | 1987-1993 | 1994-1998 | 1999-2003 | 2004-2011 | 2012-2018 
Tesechoacán 158.4 73.83 108.36 81.1 80.46 69.17 
San Juan 162.28 67.17 73.86 100.42 56.77 61.98 


En la Figura 6 se presenta la distribución espacial del potencial 


erosivo de las HRUs por subcuenca clasificado en cuatro categorías de 


erosión: ninguna a ligera, moderada, alta y muy alta (FAO, 1980). En 


TES, la mayor erosión, especialmente notable en 1986, se produjo en la 


zona alta de la subcuenca donde se concentran los bosques templados 


sobre pendientes pronunciadas, seguida de la cuenca media con erosión 


moderada a alta, mientras que la parte baja reportó volúmenes de ligeros 


a moderados en todos los escenarios. En la subcuenca SJN, por su parte, 


las zonas de mayor erosión se localizan también en los terrenos de mayor 


elevación, aunque durante los diferentes periodos se observan pequeñas 


áreas con erosión muy alta en las zonas media y baja. En esta subcuenca 


se observa una reducción del área afectada por niveles más severos de 


erosión hídrica. 
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Figura 6. Distribución del potencial erosivo por subcuenca conforme a 
la clasificación de FAO (1980). 
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En términos porcentuales (Figura 6, gráfica de áreas TES), más del 
50 % de la cuenca TES presenta erosión moderada durante todos los 
escenarios, excepto 1986 y 1998, cuando se incrementó el área con 
erosión alta. Las áreas con erosión ligera y muy alta representaron, cada 
una, menos del 20 % de cada escenario, en promedio. En la subcuenca 
SJN (Figura 6, gráfica de áreas SJN), entre 40 y 55 % del área reportó 
niveles de erosión entre 10 y 50 ton ha"! año"!, excepto en 1986; todos 


los otros niveles se mantuvieron, en su mayoría, por debajo del 20 %. 


Discusión 


La variabilidad de la precipitación y los caudales 
observados 


Para el periodo de 1986-2018, la precipitación media fue de 1 795 mm en 
TES y 2 250 mm en SJN. Ambas subcuencas registraron los menores 


niveles de precipitación en el periodo 1994-1998, con hasta 450 mm 
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menos (en SJN) que el promedio mencionado, coincidiendo con uno de 
los periodos (1997-1998) en que el fenómeno meteorológico El Niño- 
Oscilación del Sur (ENOS) tuvo mayor impacto en el territorio nacional 
(Escobar-Briones, Bonilla, Badán, Caballero, 8 Winckell, 2001). Esta 
disminución se observó en reducciones de entre 5 y 10 % del caudal 
promedio del periodo 1994-1998 en las dos subcuencas respecto a los 
demás periodos. En contraste, la precipitación del periodo 2012-2018 en 
SJN superó en más del 50 % la media de los 34 años analizados con 
máximos mensuales atípicos por encima de 500 mm en las épocas de 
lluvias de 2017 y 2018 (Conagua, 2018). 


El caudal promedio observado para los seis periodos evaluados 
osciló entre 151 y 196 m?s! en TES, y entre 208 y 251 m*s1 en SIN. 
Esta variabilidad en cada cuenca es explicada en dos escalas temporales. 
Primero, la estacionalidad de la precipitación anual se refleja en la 
amplitud de los rangos mensuales de los valores de caudales observados 
y simulados para cada periodo evaluado por subcuenca (Figura 4 y Figura 
5). Segundo, la influencia de ENOS es responsable por reducciones de 
alrededor del 10 % de la lluvia en verano e incrementos cercanos al 20 
% en la precipitación invernal en los años 1986-1987, 1991-1992, 1994, 
1997-1998 y 2002 (Magaña, Pérez, 8 Conde, 1998; Magaña, Vázquez, 
Pérez, 8i Pérez, 2003). 
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El cambio de uso de suelo y la producción de 


sedimentos 


De 1986 a 1993, TES y SIN mostraron una reducción de más del 50 % 
del volumen de erosión producida (Tabla 6), atribuible al incremento de 


la cobertura de bosques y selvas (Figura 2). 


En el periodo de 1998 se observa un aumento en el área con mayor 
severidad de erosión —de moderada (10 ton ha”! año *?) a alta (10-50 ton 
ha"? año"*)— en la zona media de la subcuenca TES (Figura 6). Esta 
agudización se relaciona principalmente con el avance de los pastizales, 
que si bien siguieron ampliándose hasta 2003 (Figura 2), ya no reflejaron 
el aumento en la erosión promedio de ese periodo. Por el contrario, en el 
escenario 2003, se observa en TES una reducción del área afectada por 
erosión alta y en una producción de sedimentos 25 % menor respecto al 
periodo anterior (Tabla 6). Esta disminución posiblemente está asociada 
con la reducción de la extensión agrícola en favor de la recuperación 
forestal, en particular en la cuenca alta, de manera semejante a lo 


encontrado por Muñoz-Villers y López-Blanco (2008) en La Antigua, 
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cuenca de la misma región hidrológica, en el mismo periodo. No obstante, 
el escenario de vegetación y uso de suelo de 2003 presentó los cambios 
más drásticos en vegetación primaria en favor de los pastizales (Figura 
2). En SIN también es marcada la expansión de los pastizales, 
especialmente en las selvas húmedas, y este cambio se tradujo en un 
incremento de la erosión en 2003 de alrededor de 35 % respecto al 


periodo anterior (Tabla 6). 


En el periodo de 2004 a 2018, TES y SJN mostraron menores 
volúmenes de erosión, posiblemente relacionados con la rápida 
recuperación típica de la vegetación tropical en las zonas alta y media 
(Figura 4), pero también se expande el área agrícola en forma de parches 


a lo largo de las dos subcuencas (Figura 2). 


Se observan ciclos de sustitución de vegetación primaria por áreas 
de uso agrícola y/o pecuario que pueden ser temporalmente abandonadas 
a la regeneración natural y luego reabiertas para el pastoreo o la 


agricultura (Alavez-Vargas et al., 2021). 


En este estudio ha sido posible relacionar los cambios en la 
producción de sedimentos con la cobertura solo cuando la cobertura 
vegetal es modificada en áreas extensas de baja pendiente, como el 
incremento de pastizales en las zonas bajas o cuando sucede en las zonas 
escarpadas de las cabeceras de cuenca, como la deforestación para uso 


agrícola. Un estudio en la cuenca La Antigua documentó, en efecto, mayor 
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producción de sedimentos en pastizales respecto a las áreas forestales 
(Martínez et a/., 2009). 


En estos casos de estudio se encontró que la topografía del 
Papaloapan está dominada por áreas con pendientes inferiores a 15 %, 
de manera que mientras las cabeceras pueden alcanzar pendientes 
fuertes a más de 2 000 msnm, más de la mitad del paisaje de cada 
subcuenca consiste en amplias planicies aluviales que funcionan como 
zonas de depósito de los sedimentos (Verstraeten et a/., 2009; Wohl et 
al., 2012). Es decir, que con independencia de su fuente, los sedimentos 
producidos no llegan a exportarse fuera de la subcuenca, pues no se 
observaron tendencias estadísticamente significativas en el volumen de 
sedimentos transportados en el cauce entre 1986 y 1998 (Alavez-Vargas 
et al., 2021). De manera semejante, el metaanálisis realizado por Guzha 
et al. (2018) en el trópico africano reveló bajas correlaciones entre la 
cobertura forestal y la escorrentía, por lo que enfatiza la necesidad de 
mejorar los programas de monitoreo de cuencas a largo plazo para 


aumentar la comprensión de las respuestas hidrológicas. 


La erosión potencial y la escasez de mediciones 
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La producción anual de sedimentos estimada con la MUSLE utiliza la 
escorrentía como factor de energía de desprendimiento y transporte de 
las partículas en sustitución del factor de erosividad de la precipitación de 
la USLE (Williams € Berndt, 1977). Por esta razón, la disponibilidad de 
datos observados para la calibración es fundamental, pues el ajuste del 
caudal implica automáticamente un ajuste a la erosión estimada y su 


omisión puede llevar a sub o sobreestimaciones significativas. 


La falta de registros históricos de mediciones de erosión se preveía 
como una restricción a los alcances de este estudio, pues limitaría la 
validación de las estimaciones. Por lo tanto, se procedió como sigue. Para 
el último periodo simulado, de 2012 a 2018, se encuentra disponible la 
cartografía de erosión de suelos de México (INEGI, 2014), elaborada 
mediante diferentes técnicas de análisis de imágenes satelitales, así como 
con información de perfiles de suelo y que presenta cuatro categorías de 
erosión basadas en los atributos de los suelos observados (Bolaños- 
González et al., 2016). Bajo los criterios de tal información, en TES, la 
mayoría de las áreas erosionadas se clasificaron de grado leve y 
moderado, presentando mayor extensión en las zonas medias y bajas, 
donde sus zonas de erosión “fuerte” (exposición del lecho rocoso en 50- 
90 % del polígono de erosión) coinciden con nuestros resultados de 


erosión alta (50-200 ton ha! año”!). 
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En SIN, los polígonos de erosión moderada se concentran en las 
zonas altas y bajas en las mismas áreas que clasificamos como de erosión 
muy alta (>200 ton ha? año**). Dicha comparación está limitada por la 
ausencia de unidades volumétricas en la cartografía referida, sin 
embargo, la cartografía de Borrelli et al. (2017) estimada a una resolución 
espacial de 250 m identifica numerosos puntos de erosión hídrica superior 
a las 20 ton ha"! año -! en nuestra área de estudio. En contraste, el estudio 
de Montes-León et al. (2011) clasifica la mayor parte de la cuenca del 
Papaloapan (y, por lo tanto, de las dos subcuencas aquí estudiadas) en 


categorías de erosión por encima de las 50 ton ha"! año”!, 


La divergencia entre las estimaciones de los estudios mencionados 
junto con el alto nivel de erosión diagnosticado en el país, referido al inicio 
de este documento, ratifican la importancia de tener información de 
campo para dimensionar mejor no solo el estado actual de los suelos y la 
cobertura vegetal de las cuencas sino también los alcances de los 
resultados de estudios como el presente. De manera similar, Tan et al. 
(2019) concluyen de su metaanálisis en el sureste asiático la urgencia de 
dirigir los esfuerzos a la generación de datos en campo de buena calidad, 


fundamentales para fines de modelación y predicción. 
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Las fuentes de incertidumbre y los alcances de las 


modelaciones 


El software SWATCUP se desarrolló como una herramienta de calibración 
para SWAT, con el objetivo de sistematizar, a través de cinco algoritmos 
posibles, los análisis de sensibilidad e incertidumbre y la parametrización 
involucrados en el proceso de calibración, dada la gran cantidad de 


parámetros incluidos en cada modelación de SWAT (Abbaspour, 2015). 


En una comparación entre los métodos de calibración de SWATCUP 
realizada en el río Ganges, India, Shivhare, Dikshit y Dwivedi (2018) 
documentaron un mejor desempeño de SUFI-2 respecto a GLUE o 
ParaSol, respaldando su amplia utilización (Abbaspour et a/., 2015; Arnold 
et al, 2012; Chotpantarat € Boonkaewwan, 2018). En esta investigación, 
la selección del grupo de parámetros óptimo para la calibración se hizo a 
partir de las recomendaciones iniciales del algoritmo de calibración y del 
análisis de sensibilidad, así como de la guía de trabajos previos 
(Abbaspour et al., 2015; Arnold et al., 2012). El grupo de parámetros 
óptimo para cada modelación incluyó lo siguiente: un parámetro de 
escorrentía superficial, tres de agua subterránea, uno de agua del suelo, 
y de 1 a 3 de evapotranspiración potencial y real (Tabla 3, Tabla 4 y Tabla 


5), mismos que han sido usados para la calibración hidrológica en otros 
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estudios recientes tanto a escala regional en Tailandia (Chotpantarat 8 
Boonkaewwan, 2018) como a escala continental en Europa (Abbaspour et 
al., 2015) con buenos resultados. Tan et al. (2019) observaron buenos 
niveles de desempeño de SWAT en la simulación de caudales en regiones 


tropicales a partir del análisis de 126 estudios del sureste asiático. 


De acuerdo con Moriasi et al. (Moriasi et a/., 2007; Moriasi, Gitau, 
Pai, € Daggupati, 2015) los valores de eficiencia de NSE > 0.50 y PBIAS 
+ 15 % pueden considerarse satisfactorios en el desempeño de las 
modelaciones de caudal mensual. Bajo estos criterios, los 12 periodos 
evaluados presentaron porcentajes de sesgo bastante buenos (entre -8.3 
y 7.5) mientras de acuerdo con el criterio de eficiencia de Nash-Sutcliffe, 
6 de 10 periodos calibrados fueron satisfactorios, siendo la excepción el 
periodo 1999-2003 en las dos subcuencas, como puede observarse en la 
Figura 4, inciso d, y Figura 5, inciso d. Tan importante como que el valor 
de la función objetivo Kling-Gupta se acerque a 1 es que el error de los 
modelos tienda a 0, es decir, que el mayor porcentaje de los datos 
observados se sitúe en el intervalo de predicción (95PPU) más estrecho 


(R-factor>0) de la mejor parametrización (P-factor>1). 


Con esta combinación de criterios, las mejores simulaciones de este 
estudio alcanzaron de forma simultánea valores de KGE y P-factor 
superiores a 0.70, y corresponden a los periodos calibrados de 1994- 
1998, 2004-2011 y 2012-2018 en TES (Figura 4, incisos c, d y e), y 1987- 
1993, 1994-1998 y 2012-2018 en SIN (Figura 5, incisos b, c, e). En un 
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rango inferior se observa el desempeño de la simulación de 1986-1993 
en TES con P-factor = 0.63 (Figura 4 inciso b), mientras que las 
simulaciones restantes presentan errores mayores al 40 %. Tal como 
documentaron Tan et al. (2019), los estudios realizados en zonas 
tropicales en condiciones de escasez de datos con frecuencia omiten la 
descripción de sus métricas de error y éstas son importantes para realizar 
una evaluación más rigurosa de los resultados obtenidos, incluso si 
evidencian, como en este estudio, anomalías importantes sobre las cuales 


prestar atención en estudios subsecuentes. 


Abbaspour et al. (Abbaspour et a/., 2007; Abbaspour et a/., 2015) 
mencionan cuatro fuentes de incertidumbre en las modelaciones: 1) la 
simplificación conceptual; 2) los procesos no simulados por el modelo; 3) 
los procesos simulados por el modelo, pero desconocidos por el modelador 
o no modelables por falta de datos, y 4) la calidad de los datos de entrada. 


En esta investigación influyeron en especial los puntos 3 y 4. 


En lo que respecta a la estimación de la erosión hídrica, se decidió 
asumir a priori la ausencia de prácticas de conservación (factor Pusie = 1 
de la MUSLE) a falta de registros históricos. En cuanto a la calidad de los 
datos, el drástico cambio en la vegetación del escenario 2003 en 
comparación con el previo (1998) y el posterior (2011), al igual que el 
bajo desempeño de las simulaciones hidrológicas de 1999-2003 en las 
dos subcuencas sugieren la necesidad de reanalizar ese escenario con un 


set de imágenes diferente. Posiblemente la irregularidad de los intervalos 
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entre los escenarios de cobertura y uso de suelo evaluados haya 


desempeñado una fuente de incertidumbre para el modelo. 


La capacidad de los modelos para reconstruir trayectorias de cambio 
hidrológico tiene como principal limitación la existencia (o no) de datos 
históricos observados para validar los escenarios generados; no obstante, 
estos continúan siendo herramientas útiles para establecer líneas base 
(Arheimer 8 Lindstróm, 2019) y, en estos casos, la comparación con 
estimaciones de otros modelos puede utilizarse como estrategia de 
validación de los resultados (Borrelli et a/., 2017) y de delimitación de los 


alcances. 


Conclusiones 


Se estudió la respuesta de la erosión potencial en dos subcuencas 
hidrológicas tropicales de mediana escala respecto a los cambios de uso 
de suelo. La herramienta para llevar a cabo dicha cuantificación fue el 
modelo hidrológico semidistribuido Soil and Water Assessment Tool 
(SWAT). 
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De los 12 periodos simulados, dos fueron invalidados por anomalías 
en el escenario de uso de suelo. Basados en la evaluación combinada de 
tres estadísticos de eficiencia y dos de incertidumbre, de los 10 periodos 
calibrados restantes, seis mostraron buenos ajustes en la simulación 
hidrológica con KGE > 0.70 y P-factor > 0.70. De acuerdo con las 
estimaciones realizadas, se observó una reducción de los volúmenes de 


erosión hídrica a medio plazo. 


La calidad de los resultados se evaluó a través del uso combinado 
de tres estadísticos de eficiencia y dos medidas de la incertidumbre. Al 
usar modelos hidrológicos para la investigación y el soporte a la toma de 
decisiones de política pública se recomienda a los profesionales y 
académicos la selección de una combinación de estadísticos equivalente 
a la empleada en este estudio, con el fin de verificar de manera conjunta 


y con mayor rigor la confiabilidad de las modelaciones. 


La comparación de los volúmenes de erosión potencial estimada 
respecto a los resultados obtenidos por otros estudios en la misma área 
de estudio, pero con métodos y escalas diferentes, permite validar los 


resultados y acotar los alcances. 


Asimismo, la heterogeneidad observada entre los estudios 
revisados y el alto nivel de erosión diagnosticado en el país permite 
argumentar en favor de la necesidad de generar estándares y criterios 
para la producción de datos de buena calidad para la investigación y 


gestión de la erosión en el territorio nacional. 
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En países como México, con monitoreo limitado, la estimación de 
las variables de interés a partir de datos físicos, climáticos y humanos, así 
como la disponibilidad de bases de datos de remuestreo, ha popularizado 
el uso de SWAT. Las funcionalidades de esta herramienta son valiosas 
para prospectar cuencas no aforadas o cronosecuencias de interés, sin 
embargo sus aplicaciones para fines de investigación o gestión deben 
estar sustentadas tanto con métricas de error como con estadísticos de 
eficiencia que dimensionen con claridad los alcances, limitaciones y 


anomalías de los resultados obtenidos. 


Los modelos de base física son herramientas de utilidad en las 
regiones tropicales con escasez de datos históricos para estimar el 
comportamiento de las dinámicas hidrológicas y edafológicas, pero dados 
los niveles de incertidumbre encontrados en este trabajo, se considera 
que no pueden suplir de manera permanente las funciones de una red de 


monitoreo de erosión a escala de cuenca. 
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